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On  singular  values  of  Hankel  operators  of  finite  rank 

W.  B.  Graggt  and  L.  Rcichel  | 


Abstract 

Let  Hbtz  Hankel  operator  defined  by  its  symbol  p  =  k/x  where  x  is  a  monic  polynomial  of  degree 
n  and  IT  is  a  polynomial  of  degree  less  than  n.  Then  H  has  rank  n.  We  derive  a  generalized  Takagi 
singular  value  problem  defined  by  two  n  x  n  matrices,  such  that  its  n  generalized  Takagi  singular 
values  arc  the  positive  singular  values  of  H.  If  p  is  real,  then  the  gencr^ized  Takagi  singular  value 
problem  reduces  to  a  generalized  symmetric  eigenvalue  problem.  The  computations  can  be  carried 
out  so  that  the  Lanczos  method  applied  to  the  latter  problem  requires  only  0(n  log  n)  arithmetic 
operations  for  each  iteration.  If  -k  and  x  given  in  power  form,  then  the  elements  of  all  n  x  n 
matrices  required  can  be  determined  in  O(n’)  arithmetic  operations. 
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1.  Introduction 


Let  H  =  (»7y+fc]^=o  be  a  be  a  Hankel  operator  defined  by  its  rational  symbol  p  —  jt/x,  where 

n  — 1  n 

j=0  ]=o 


We  assume  that  tt  and  x  have  no  common  zeros.  The  elements  »?y  of  H  are  then  given  by 


^w  =  i|  =  EM-'-.  (1.2) 

In  order  to  simplify  our  presentation,  we  assume  that  the  zeros  of  x  are  distinct.  How  our 

formulas  need  to  be  modified  in  order  to  remove  this  assumption  is  discussed  in  Remark  1.1  below. 
Hence  p  has  a  partizd  fraction  decomposition 

KA)=:Ex^.  (1-3) 

*=1  * 

Expansion  of  the  right  hand  side  of  (1.3)  into  a  geometric  series,  and  comparison  with  (1.2),  yields 

n 

=  (1-4) 

fc=l 


We  now  express  (1.4)  in  matrix  form.  Let 


A  :  =  diag[ai,a2,...  a„]  c 

(1.5) 

A:=  diaslAi.Aj,... 

(1.6) 

and  introduce  the  Vandermonde  matrix 

Vo  :  =  lAi+xl7.*io  e 

(1.7) 

Define 

V  := 

(1.8) 

where 

V,  :=  VoA^'‘,  j  >  1. 

(1.9) 

Then  (1.4)  can  be  written  as 

H  =  VAV'^. 

(1.10) 

Let  denote  the  vector  space  C°°  equipped  with  the  Euclidean  norm. 


2 


Proposition  1.1.  H  —*  P  bounded  |Afcl  <  1  for  1  <  k  <  n. 

Proof.  The  proposition  holds  independent  of  the  multiplicity  of  the  A*.  In  the  present  proof  we 
assume  that  the  Ajt  are  distinct.  The  proof  for  confluent  Afc  is  commented  on  in  Remark  1.1. 

Let  cj  =  [ejI^o  ^  vector  with  eq  =  1-  Then 

^  ~  G  »7y  0  as  j  oo  => 

|Afc|  <  1  for  1  <  A:  <  n, 
where  the  last  implication  follows  from  (1.4). 

Conversely,  assume  that  lAfc|  <  1  for  1  <  A:  <  n.  Then  by  (1.8)  -  (110)  we  obtain 

iiffib  <  uh\\v\'&  <  iiAii,iiv„iiiiif;(A''Ar'ig 

j=o 

=  l|A|bl|V'olllll(/-(A''An-M|3.  ■ 

We  assume  henceforth  that  |Afc|  <  1  for  1  <  A:  <  n.  Introduce 

U:^VV^\  (1.11) 

Ho~VoAV^.  (1.12) 

Then  Hq  has  rank  n.  We  note  ,  by  comparing  (1.12)  with  (1.10),  that  Ho  is  the  leading  principal 
n  X  n  submatrix  of  H.  From  (l.lO)  -  (112)  it  follows  that 

H^UHoV'^.  (1.13) 

The  leading  n  x  n  submatrix  of  U  is  the  n  x  n  identity  matrix.  U  therefore  is  of  rank  n  and 
can  be  factored 

U  =  QR,  i?  e  (T"''", 

where  Q^Q  =  In  and  i?  is  a  nonsingular  right  triangular  matrix.  We  obtain 

(7+{H)  =  a+{QRHoR^Q'^)  =  o{RHoR^),  (1.14) 

where  a  denotes  the  set  of  singular  values  and  denotes  the  subset  of  the  positive  ones. 

The  n  X  n  matrix  RHqRT  is  complex  symmetric.  Takagi  [Tal],  lTa2]  showed  the  existence  of  a 
complex  syrrunetric  singular  value  decomposition 

RHqR'^  =W'EW'^ ,  WeC"’'",  E  =  diag\au(^2,.:  <Tn],  (115) 

where  W^W  =  /„  and  <7y  >  0  are  the  singular  values  of  RHqRT  .  In  Section  2  we  present  an 

elementary  proof  of  the  existence  of  this  decomposition.  Let  W  =  [lui ,  ii;2,  •••  Then 

(1.15)  can  be  written  Jis  the  Takagi  singular  value  problem 

RHoR^Wj  —  Wj(Tj,  Wk  =  1  <  J,  A:  <  n,  (116) 
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where  the  bar  denotes  complex  conjugation  and  is  Kronecker’s  6  function.  The  problems  (1.15) 
-  (1.16)  could  be  solved  by  the  algorithm  described  in  [BGG],  but  this  would  require  RHqBT  to  be 
explicitly  computed.  In  order  to  avoid  these  matrix  multiplications  we  let  Vj  :=  R^Wj  and  obtain 
from  (1.16)  the  generalized  Takagi  singular  value  problem 

HqVj  =  {R” R)~^Vj(7j,  v^{R“ R)~^Vk  =  Sjk,  1  <j,k  <  n.  (1.17) 

The  solution  of  (117)  requires  {R^ R)~^  to  be  known.  In  Section  3  we  show  that 

{R^R)-^  =  I  -  BoB^,  (1.18) 

where  Bq  e  is  a  triangular  Toeplitz  matrix.  The  elements  of  Bq  and  Ho  can  be  determined 

from  the  coefRcients  of  jt  and  x  0(n  log  n)  arithmetic  operations  by  the  fast  Fourier  transform 
(FFT)  method.  This  is  demonstrated  in  Section  4.  Section  5  shows  that 

R^R  =  TiMoT",  Ti.Afo  e  (1.19) 

where  Ti  and  Mo  are  Toeplitz  matrices,  and  describes  a  numerical  scheme  for  the  computation 
of  this  factorization  from  (1.16)  in  O(n^)  arithmetic  operations.  We  also  present  a  Hermitian 
factorization  of  R^  R  into  n  x  n  triangular  matrices. 

The  factorization  (1.19)  may  be  of  interest  for  the  numerical  solution  of  (1.17).  Assume  that 
the  coefficients  of  ir  and  x  are  real  valued.  Then  Hq,  {R^ R)~^  e  and  (117)  reduces  to 

a  generalized  symmetric  eigenvalue  problem.  The  Lanczos  method  ([Pa,  Section  15.11),  [ER]) 
would  appear  suitable  for  solving  this  eigenproblem  for  the  following  reason.  Let  C  e  be  a 

Hankel  or  Toeplitz  matrix  and  let  v  e  C"  be  arbitrary.  It  is  well  known  that  Cv  can  be  computed 
in  0(n  log  n)  arithmetic  operations  using  FFTs.  Hence  Hqv,  {R^ R)~^v  and  [R^  R)v  can  be 
computed  in  0(n  log  n)  arithmetic  operations,  where  we  use  (118)  -  (1.19).  Each  iteration  of  the 
Lanczos  algorithm  given  in  [Pa,  p.324|  therefore  requires  only  0(n  log  n)  arithmetic  operations. 

The  computation  of  singular  values  of  H  is  important  in  Hankel  norm  approximation  problems  of 
systems  theory,  such  as  the  model  reduction  problem  [Gl).  The  approximation  of  functions  by  the 
Caratheodory  -  Fejer  method  yields  another  application  [Gu],  [IV]. 

Other  methods  for  reducing  the  singular  value  problem  for  H  to  a  finite  dimensional  one  have  been 
described  by  Kung  and  Gutknecht  [Gu]  and  Young  [Yo].  These  methods,  however,  do  not  preserve 
symmetry.  Moreover,  Young’s  approach  requires  generally  O(n^)  arithmetic,  operations  to  compute 
the  matrices  required. 


Remark  1.1.  Formulas  (1.3)  -  (1.8)  and  the  proof  of  Proposition  1.1  require  distinct  Afc.  This 
restriction  can  be  removed.  Assume  first  that  Aj  =  Aj  =  •••  =  ■^n.  Then  (1.3)  -  (1.4)  have  to  be 
replaced  by 

n 


pW  =:  E 

fc=l 


Qk 

(A -A,)*’ 


(1.3') 


= 


E 


Qk 


k  1  ^ 

«=!  j=o 


>’1' 


(1.4') 


In  (1.5)  A  has  to  be  substituted  by  the  upper  triangular  Hankel  matrix 


A  —  [aj+fc+i]"fc=o  f  ^ 


Op  :=  0,  p  >  n. 
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The  matrix  A  in  (1.6)  has  to  be  replaced  by  the  Jordan  matrix  with  all  diagonal  elements  equal  to 
Ai  and  all  superdiagonal  elements  equal  to  one.  The  matrix  Vq  in  (l.V)  need  be  replaced  by  the 
onfluent  Vandermonde  matrix.  For  instance,  we  obtain  for  n  =  3 


With  A,  A  and  Vq  modified  as  described,  we  define  Vy  and  V  by  (1.8)  -  (1.9),  U  by  (1.11)  and 
Hq  by  (1.12).  Then  (1.10)  and  (1.13)  hold  and  Hq  is  the  leading  principal  n  x  n  submatrix  of  H. 
Also  (1.14)  -  (1.19)  remain  valid.  Proposition  1.1  can  be  shown  by  replacing  (1.4)  by  (1.4'),  and 
by  bounding  the  sum 

II  £(*"*)"' III 

j=o 

where  A  now  is  a  Jordan  matrix.  This  sum  is  bounded  if  jAj  |  <  1,  and  the  proposition  remains 
valid. 

In  general,  when  the  Xk  are  of  arbitrary  multiplicity,  A  in  (1.5)  has  to  be  replaced  by  a  block 
diagonal  matrix,  where  each  block  is  an  upper  triangular  Hankel  matrix.  The  blocks  are  of  the 
same  sizes  as  the  multiplicities  of  the  A*,  and  the  number  of  blocks  equals  the  number  of  distinct  Afc. 
A  in  (1.6)  is  replaced  by  a  Jordan  matrix  with  Jordan  boxes  of  the  same  sizes  as  the  multiplicities 
of  the  Afc,  and  the  number  of  boxes  equal  to  the  number  of  distinct  Afc.  Vq  in  (1.7)  is  replaced  by 
an  appropriate  confluent  Vandermonde  matrix.  With  these  changes  (1.10)  -  (1.19)  are  valid,  and 
so  is  Proposition  1.1.  We  omit  the  details  since  the  numerical  computations  are  independent  of  the 
multiplicity  of  the  Afc.  « 
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2.  The  Symmetric  Singular  Value  Decomposition 


In  this  section  we  present  an  elementary  proof  of  Takagi’s  theorem,  i.e.  we  show  the  existence  of 
a  symmetric  singular  value  decomposition  of  a  complex  symmetric  matrix.  Let  C  =  c 
and  define  A,  B  e  by  C  :=  A  +  iB,  i  :=  v/— 1.  Then  A  =  A^and  B  =  B^ ,  so  the  matrix 


is  real  and  symmetric.  Let  be  the  positive  eigenvalues  of  C  and  form 

E  :=  diag\ax,a2y  -,Or\- 
Let 


with 

and 

Write  (2.1)  as 


A 

B 

U' 

u' 

B 

-A 

V 

V 

{ 


U,Ve 


U'^U  +  V^V  =  Ir. 


AU  +  BV  =  UE 
BU  -  AV  =  VE 


and  note  that  (2.2)  also  can  be  written  as 

{ 


i.e. 


with 


AV  +  B{-U)  =  V(-E) 

BV  -  A{-U)  =  (-t/)(-E), 


(-S) 


V'^V  +  {-Uf{-U)  =  /,. 


A 

B 

V' 

v' 

B 

-A 

-U 

-U 

(2.1) 


(2.2) 


(2.3) 


Hence  C  has  at  least  r  negative  eigenvalues.  We  could  also  have  let  ay  be  the  negative  eigenvalues  of 
C  and  then  (2.3)  would  have  given  us  positive  ones.  We  therefore  may  assume  that  ±ai ,  ±a2, ...,  ±ar 
are  ail  the  nonzero  eigenvalues  of  C. 


Since  eigenvectors  associated  with  distinct  eigenvalues  of  a 
we  have 


0  = 


=  V^U 


real  symmetric  matrix  are  orthogonal, 
-U'^V. 


The  spectral  resolution  of  C  is  thus 


A  B' 

U  V 

E 

B  -A 

V  -U 

-E 

6 


which  yields 


r  A  =  C/EU^  -  VEV'^ 
\b  =  VEU^  +  UEV^. 


Therefore 


where 


Moreover 


C  =  A  +  iB  =  UEU"^  -  VEV^  +  i{VEU'^  +  UEV'^) 

r 

=  {U  +  iV)E{U‘^  +  »V^)  =  WEW'^  = 

k  =  l 

U  +  iV  =:  W  =  (T". 


W^W  =  {U'^  -  iV^)(C/  +  iV)  =  {U'^U  +  V'^V)  +  i{U^V  -V’^U)  =  Ir. 
If  r  <  n  then  one  may  replace  E  by 


and  W  by 


where  tWr  +  i , tUnC  C" 


Eo  ;=  diag[(Ti,<T2,...,(rr,0,  •••,0]  e 

Wo  :=  [u;:,u;2,-.,«^r,ti;r+i,  •••u'n]  f  C"’'", 
are  chosen  so  that  W^Wq  =  In  ■ 
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3.  A  SLiiple  Expression  for  [R^  R)  ^ 

In  this  section  we  derive  (1.18).  Introduce  the  Frobenius  matrix 

F:= 

where 

cy  :  =  <  R",  2  <  i  <  n,  (3.1) 

/  :  =  [xo.Xi.  ".Xn-i)’'<  4^"- 

Then  F  is  the  companion  matrix  of  x  and 

F^Vo  =  VoA.  (3.2) 

Throughout  this  section  Vq  and  A  are  defined  by  (1.6)  -  (1.7)  if  the  A*  are  distinct.  For  confluent 
Afc  we  modify  Vq  and  A  according  to  Remark  1.1.  The  following  lemma  shows  that 

G  (3.3) 

satisfies  a  Stein  equation.  This  will  enable  us  to  obtain  a  simple  expression  for  by  an  application 
of  the  Sherman-Morrison- Woodbury  formula. 


Lemma  3.1.  G  is  the  unique  solution  of  the  Stein  equation 


X  -  =  /„, 


Xe 


Proof.  By  (1.8),  (1.9)  and  (1.11)  we  obtain 


R^R  =  U^U  =  ^  Vo“"(A"*)"Vo"VbA’’''Vo“\ 

k=0 


and  (3.2)  yields  now 

OO 

(7=5^  F'**(F"*)". 

fc=0 


(3.4) 


(3.5) 


(3.6) 


The  series  in  (3.5)  -  (3.6)  converge  because  lAfc|  <  1  for  all  k.  Substitution  of  (3.6)  into  (3.4) 
shows  that  G  solves  (3.4).  The  unicity  follows  from  |Afc|  <  1  for  all  k.  The  latter  can  be  seen  by  a 
similarity  transform  of  F"  to  Schur  triangular  form,  a 


Introduce  the  cyclic  downshift  operator  in 


where 


Let 


h.ea-  .en.e,] 

:=  . fR"". 

t  ~  Ixo.Xi-  ■.Xn,0,0,...,0]^c 


(3.7) 
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and  define  the  Toeplitz  matrix  T  of  parallelogram  form 

T  :=  e  (3.8) 

Let  To  be  the  leading  n  x  n  submatrix  of  T,  and  let  T\  be  the  trailing  n  X  n  submatrix  of  T.  Then 
To  is  a  left  triangular  Toeplitz  matrix,  and  Ti  is  a  unit  right  triangular  Toeplitz  matrix. 

Example  3.1.  Let  n  =  3.  Then 
Xc 

Xi  Xo 
X2  X2  Xo 
X3  X2  Xl 
X3  X2 
.X3 

where  we  note  that  ;)f3  =  1.  • 

Lemma  3.2.  Let  To  and  Ti  be  defined  as  above.  Then 

To"  To  +  T"  Ti  =  ToTo"  +  Tj  T" .  (3.9) 

Proof.  Let  N  ;=  —  T"To  +  T"Ti.  We  first  show  that  N  is  &  Toeplitz  matrix.  Let  Cj  be 

defined  by  (3.1).  Then  by  (3.8)  we  have  for  1  <  j,k  <  n, 

ej  Nck  =  eyT"T  e*  =  E’^'H, 

where  we  have  used  that  E^  —  E~^ .  We  next  define  the  reversal  matrix 

J:= 

Toeplitz  matrices  are  counter  symmetric,  i.e.  N  =  JN^  J.  Using  that  N  is  counter  symmetric  and 
Hermitian  yields 

To"  To  +  Ti"Ti  =  ^■  =  JN'^J  =  JNJ  =  J{T^%  + 

=  JT^  J  •  J%J  +  JTJ  J  •  JT^J  =  ToTo"  +  TiTj" .  « 

The  next  lemma  presents  a  Gaussian  factorization  of  F"  in  terms  of  To  and  T\.  This  will  be  used 
together  with  Lemma  3.1  to  express  G~^  in  terms  of  Tq  and  Tj. 

Lemma  3.3. 

F"  =  -ToTf^  (3.10) 

Proof.  We  first  show  that 

r  Vo 

(To^.T^  =0.  (3.11) 

VoA- 

Let  ey  be  defined  by  (3.7)  and  assume  for  the  moment  that  the  are  distinct.  Then 

r  Vo 

^][TLtT\  «*  =  x(Afc)Ar'  (3.12) 

Vo  A” 


Xo  X3  X2  Xl 

To  =  Xl  Xo  ,  Ti  =  X3  X2  , 

.X2  Xl  XoJ  [  X3. 


9 


and  the  right  hand  side  vanishes  for  I  <  j,k  <  n.  If  the  A*  are  confluent,  then  the  right  hand  side 
expression  of  (3.12)  contains  derivatives  of  x{A)  evaluated  at  Xk-  The  right  hand  side  of  (3.12), 
however,  still  vanishes  and  (3.11)  holds. 

We  now  write  (3.11)  as 

To^Vb  +  TfVoA"  =  0 
and  apply  (3.2).  This  shows  (3.10).  ■ 

We  are  now  in  a  position  to  show  (1.18).  By  (3.4)  G  satisfles 

G  =  I  +  F''G 

and  an  application  of  the  Sherman-Morrison- Woodbury  formula  yields 

G-i  =  (/  +  F"-^)~^  =  I  -  F"(G"^  +  (3.13) 


We  now  determine  an  expression  for 

Y:=I-G-\  (3.14) 

Substitute  Y  zuid  (3.10)  into  (3.13)  to  obtain 


Y  =  To{T"To  +  T{^Ti  -  T"YTi)-^T^.  (3.15) 

In  order  to  determine  a  simple  expression  for  Y  from  (3.15)  we  need  the  following  observation, 
which  is  also  central  to  Section  4.  To  and  are  both  left  triangular  n  x  n  Toeplitz  matrices. 
Multiplication  of  To  with  can  be  identified  with  polynomial  multiplication,  see  [Hel,  Section 
1.3]  and  Section  4.  Since  multiplication  of  polynomials  commutes,  we  obtain 

ToTf"  =  Tf"To.  (3.16) 


From  the  correspondence  between  polynomials  and  left  triangular  Toeplitz  matrices  it  also  follows 
that  is  a  left  triangular  Toeplitz  matrix. 

Lemma  3.4.  Equation  (3.15)  has  the  unique  solution 

Y  =  TqTq" Tj- ^  =  ToTf  "  Tf  *  To"  .  (3.17) 

Proof.  Unicity  follows  from  (3.14)  and  that  (3.4)  has  a  unique  solution.  From  (3.16)  we  obtain 

Tf^ToTo^Tf^  =ToTi-"Tf^To".  (3.18) 

Now  substitute 

K  =  Tf  "  TqT"  Tf  ^ 

into  (3.15).  We  obtain 

Tr"ToTo"Ti-^  =  To(To"To  +  T"Ti  -  ToTo")-'To".  (3.19) 
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An  application  of  (3.9)  reduces  (3.19)  to  (3.18).  The  latter  has  already  been  shown  to  be  valid. 
Therefore  (3.17)  solves  (3.15).  ■ 


Let 

Then  Bq  is  a  left  trizuigular  n 


Bo  :=  ToTf’'  =  Tr^To. 
n  Toeplitz  matrix.  By  (3.14)  and  (3.17) 

G-^  =  /  -Bo^o  =I-B^Bo. 

{R^ R)-^  =  I  -  BoB^ . 


(3.20) 


From  (3.3)  it  now  follows  that 


(3.21) 


4.  Computation  of  Hq  and  Bq 

We  summarize  some  results  in  (He  1,  Section  1.3]  and  [He  2,  Section  13.9]  in  order  to  show  that  the 
elements  of  Hq  and  Bq  can  be  computed  in  0(n  log  n)  arithmetic  operations  from  the  coefficients 
Xj  of  X  and  jry  of  jt,  see  (11).  To  a  polynomial  or  power  series 

fW  -E  &  a'  +  o(a") 

3=0 

we  associate  the  left  trizuigular  n  x  n  Toeplitz  matrix 

Z  J  <  0> 

and  we  write  f  — »  Z.  If  ^(A)  is  a  polynomial  and  X  a  left  triangular  n  x  n  Toeplitz  matrix  such 
that  ^  X,  then  it  is  easily  seen  that  — »  ZX.  In  particular,  ZX  is  a  left  triangular  n  x  n 
Toeplitz  matrix.  From  and  — ►  XZ  is  follows  that  ZX  =  XZ. 

Assume  that  fo  0  and  let  1/f  — »  Z' .  Then  1/f  •  f  — +  /,  Z'Z  and  ZZ' .  We  obtain  Z'  =  Z~^  and 
therefore  Z~^  is  a  left  triangular  Toeplitz  matrix. 

Example  4.1.  We  have  x  '^o-  Let 


i;(A)  :=  A"  x(l/A)  =  E  An-,A>-  (4.1) 

3=0 

Then  x  and  the  Blaschke  product 

4  ^ToT{^  =  Bo.  (4.2) 

X  ■ 

Now  let  ^(A)  and  f(A)  be  jtfbitrary  polynomials  such  that  f(0)  ^  0.  Henrici  [He2,  Theorem  13.9e  ] 
shows  that  the  first  n  coefficients  in  the  MacLaurin  expansion  of  ^(A)/f(A)  can  be  computed  in 
0(n  log  n)  multiplications.  The  proof  uses  FFT.  It  is  easily  seen  that  the  number  of  additions  also 
is  0(n  log  n). 

From  Xr»  =  1  and  (4.1)  we  obtain  x(0)  0.  Hence,  the  first  n  terms  in  the  MacLaurin  expansion 

of  x/x  can  be  computed  in  0(n  log  n)  arithmetic  operations.  By  (4.2)  therefore  ToTf^  =  Bq  can 
be  computed  in  0{n  log  n)  arithmetic  operations. 


Because  A"x(l/^)  ^  0  for  A  =  0,  we  can  compute  the  first  n  terms  in  the  MacLaurin  expansion  of 


A"x(l/A) 

A'‘X(1/A) 


=  E-);V^'+0{A") 

3=0 


in  0(n  log  n)  arithmetic  operations.  This  shows  that  Hq  can  be  computed  in  0(n  log  n)  arithmetic 
operations. 
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5.  A  Factorization  of  R 

It  follows  from  (3.3)  and  (3.20)  -  (3.21)  that 

G-^  =  {R^R)~^  =  1-  Bo  =  /  -  Tf  ^^oTo"T^^  (4.1) 

and  therefore 

T^G-^Ti  =  T^Ti  -  ToT^  =:  MqK  (4.2) 


The  expression  defining  Mq^  is  a  Gohberg-Semencul  formula  for  the  inverse  of  an  n  x  n  Toeplitz 
matrix,  see,  e.g.,  (lo.  Theorem  18.2,  p.  152].  We  denote  this  Toeplitz  matrix  by  Mq.  From  the  left 
hand  expression  of  (4.2)  and  the  nonsingularity  of  Ti  and  R  it  follows  that  Mo  is  Hermitian  and 
positive  definite.  The  desired  factorization  of  R^ R  is 


R^R  =  Ti  Mo  T". 


We  will  now  show  how  Mq  can  be  computed.  The  computation  involves  running  the  Levinson 
algorithm  backwards. 

Consider  the  related  Gohberg-Semencul  formula,  see,  e.g.,  [lo.  Theorem  18.1,  p.  148]  or  [AG], 


- 

- 

H 

■ 

- 

Xn  Xn-1 

Xo 

Xn  Xn-1  ■  • 

xo 

Xn— 1 

X 

X 

Xn-1 

Xn 

Xn 

(4.3) 


■  0 

•  0 

Xo 

Xo 

Xl 

\  • 

-Xn-l  •••  Xl  Xo  0. 

-Xn-l  Xl  Xo  0. 

where  the  four  triangular  Toeplitz  matrices  define  the  inverse  of  an  (n  -|-  1)  x  (n  1)  Hermitian 
Toeplitz  matrix.  Denote  this  Toeplitz  matrix  by  Mi.  Then  Mo  is  the  leading  principal  n  x  n 
submatrix  of  Mi,  see  [lo.  Theorems  18.1  -  18.2). 

Let  Ri  :=  [/>yfc]”fc=o  ^  be  the  unit  right  triangular  matrix,  and  let 

Di  :=  diag\6o,Si,...,Sn]  be  the  diagonal  matrix  such  that 

R^MiRi^Di.  (4.4) 
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Given  Afi  =  (A*y-fc]",fc=oi  matrices  Ri  and  Di  can  be  computed  by  the  Levinson  algorithm,  and 
by  comparing  Ri  with  (4.3)  one  finds  that 

Pin  =  Xj.  0  <  y  <  n  and  =  Xn, 

see,  e.g.,  (AG).  We  now  apply  the  recursion  formula  in  Levinson’s  algorithm  backwards  in  order 
to  determine  Ri  and  Di  from  the  last  column  of  Ri  and  Then  the  recursion  formula  is  used 
forwards  to  determine  Mq.  We  will  also  obtain  a  Hermitian  factorization  of  R^  R  into  triangular 
matrices. 


Backward  Levinson  algorithm 

input:  [Pjn]y=o>^n  •  outpi  l:  Ri,Di,  Schur  paran -iers  {7y}"=i  of  Mq] 
for  A:  :=  n,n  -  l,n  -  2, ...,  1  do 

Ik  ■=  Pok\  Pk-i,k-i  :=  1; 
for  y  :=  1,2,  ...,integ<  r  part(5)  do 
1  solve  for  pj^i^k-i  and  pk-i-j,k-i  the  linear  system  of  equations 


1 

7fc 

P]-l.k-l 

Pj.k 

1  J 

.Pk-l-],k-l  . 

■  Pk—j,k  . 

:=  (W(l-|7.|))/(l  +  l7.|); 

Levinson  recursion  for  computing  Mq  =  [fij-k]^kio 
input:  i2i,I?i,{7j}"=i:  output: 

Mo  :=  ^o;  Ml  “  “^7ii 
for  k  :=  1,2,  ...,n  -  1  do 

Mfc+l  “  ~^klk+i  ~  ^i-1  MjPy-i,fei 


Hence  Mq,  Ri,  and  Di  are  computed  in  0{n^)  arithmetic  operations  from  the  coefficients  of  x-  Let 
Rq  and  Dq  denote  the  n  x  n  leading  principal  submatrices  of  i?i  and  Di  respectively.  Similarly  to 
(4.4)  we  have 

R^MqRq  =  Do.  (4.5) 

1  /2 

Because  Mq  is  positive  definite,  so  is  Dq.  Dq'  can  therefore  easily  be  computed.  Wp  obtain  from 
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(4.1)  -  (4.2)  and  (4.5),  with  R  := 

{RT^fiM'f).  (4.6) 

The  right  hand  side  of  (4.6)  is  a  Hermitian  factorization  into  triangulu  matrices.  It  ceui  be  computed 
in  0(  n^'f  arithmetic  operations  from  the  coefficients  of  x- 
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